package scu.maqiang.nla;

import scu.maqiang.fes.*;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;
import java.util.Arrays;

public class NonlinearPossion2D {
    public static void main(String[] args) {
        //NewtonIteration_HeatTransfer();
        //DirectIteration_HeatTransfer();
        //DirectIteration_example662();
        //NewtonIteration_example662();
        //DirectIteration_membrane();
        //NewtonIteration_Membrane();
        //DirectIteration_HeatTransfer_2DQ4();
        //NewtonIteration_HeatTransfer_2DQ4();
        //DirectIteration_example662_2DQ4();
        //NewtonIteration_example662_2DQ4();
        //DirectIteration_membrane_2DQ4();
        //DirectIteration_membrane_2DQ9();
        //NewtonIteration_Membrane_2DQ4();
        NewtonIteration_Membrane_2DQ9();
    }

    public static void DirectIteration_HeatTransfer() {
        Mesh2T3 mesh = new Mesh2T3().square2D(8, 4, T3Type.Left);
        mesh.scale(0.18, 0.1);
        mesh.displayNodes();
        mesh.displayElements();
        mesh.displayBoundarys();

        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 0.2;
        double k1 = 2e-3;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] k0k1Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0k1Vec, i-> k0 * k1);
        double[][] coef = {k0Vec, k0k1Vec, zeroVec, zeroVec,
                           k0Vec, k0k1Vec, zeroVec, zeroVec, x};
        while(k < maxIteration) {
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_1, A);
            fs.applyBC_MBN(A, RHS, Direct.All, 500, 4);
            fs.applyBC_MBN(A, RHS, Direct.All, 300, 2);
            solver.setMatrix(A);
            solver.PCGSSOR(RHS, dx, 1.5, 1);
            double[] diff = MVO.add(dx, -1.0, x);
            if(MVO.L2Norm(diff) <= 1.0e-8 * MVO.L2Norm(x)) {
                break;
            }
            Arrays.setAll(x, i -> dx[i]);

            A.clear();
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("FinalSol.dat", x);
    }

    public static void DirectIteration_HeatTransfer_2DQ4() {
        Mesh2Q4 mesh = new Mesh2Q4().square2D(8, 4);
        mesh.scale(0.18, 0.1);
        mesh.displayNodes();
        mesh.displayElements();
        mesh.displayBoundarys();

        FES2Q41 fs = new FES2Q41(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 0.2;
        double k1 = 2e-3;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] k0k1Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0k1Vec, i-> k0 * k1);
        double[][] coef = {k0Vec, k0k1Vec, zeroVec, zeroVec,
                k0Vec, k0k1Vec, zeroVec, zeroVec, x};
        while(k < maxIteration) {
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_1, A);
            fs.applyBC_MBN(A, RHS, Direct.All, 500, 4);
            fs.applyBC_MBN(A, RHS, Direct.All, 300, 2);
            solver.setMatrix(A);
            solver.PCGSSOR(RHS, dx, 1.5, 1);
            double[] diff = MVO.add(dx, -1.0, x);
            if(MVO.L2Norm(diff) <= 1.0e-8 * MVO.L2Norm(x)) {
                break;
            }
            Arrays.setAll(x, i -> dx[i]);

            A.clear();
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("FinalSol2DQ4.dat", x);
    }

    public static void DirectIteration_example662() {
        Mesh2T3 mesh = new Mesh2T3().square2D(8, 4, T3Type.Left);
        mesh.scale(0.2, 0.1);
        mesh.displayNodes();
        mesh.displayElements();
        mesh.displayBoundarys();

        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 0.2;
        double k1 = 100.0;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] k0k1Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0k1Vec, i-> k0 * k1);
        double[][] coef = {k0Vec, k0k1Vec, zeroVec, zeroVec,
                k0Vec, k0k1Vec, zeroVec, zeroVec, x};
        ScalarFunc topBoundary = (xy, label, param) -> 500 * (1 - 10 * xy[0] * xy[0]);
        while(k < maxIteration) {
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_1, A);
            fs.applyBC_MBN(A, RHS, Direct.All, 500, 4);
            fs.applyBC_MBN(A, RHS, Direct.All, 300, 2);
            fs.applyBC_MBN(A, Direct.All, 3);
            fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[]{topBoundary}, null, 3);

            solver.setMatrix(A);
            solver.PCGSSOR(RHS, dx, 1.5, 1);
            double[] diff = MVO.add(dx, -1.0, x);
            if(MVO.L2Norm(diff) <= 1.0e-8 * MVO.L2Norm(x)) {
                break;
            }
            Arrays.setAll(x, i -> dx[i]);

            A.clear();
            k++;
        }
        //System.out.println(Arrays.toString(x));
        mesh.toTecplot("SolExample662.dat", x);
        A.clear();
        fs.assembleStiff(new double[]{k0}, BVPType.COMMON, A);
        Arrays.fill(RHS, 0.0);
        fs.applyBC_MBN(A, RHS, Direct.All, 500, 4);
        fs.applyBC_MBN(A, RHS, Direct.All, 300, 2);
        fs.applyBC_MBN(A, Direct.All, 3);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[]{topBoundary}, null, 3);
        solver.setMatrix(A);
        solver.PCGSSOR(RHS, x, 1.5, 1);
        mesh.toTecplot("LinearSolExample662.dat", x);
    }

    public static void DirectIteration_example662_2DQ4() {
        Mesh2Q4 mesh = new Mesh2Q4().square2D(8, 4);
        mesh.scale(0.2, 0.1);
        mesh.displayNodes();
        mesh.displayElements();
        mesh.displayBoundarys();

        FES2Q41 fs = new FES2Q41(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 0.2;
        double k1 = 100.0;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] k0k1Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0k1Vec, i-> k0 * k1);
        double[][] coef = {k0Vec, k0k1Vec, zeroVec, zeroVec,
                k0Vec, k0k1Vec, zeroVec, zeroVec, x};
        ScalarFunc topBoundary = (xy, label, param) -> 500 * (1 - 10 * xy[0] * xy[0]);
        while(k < maxIteration) {
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_1, A);
            fs.applyBC_MBN(A, RHS, Direct.All, 500, 4);
            fs.applyBC_MBN(A, RHS, Direct.All, 300, 2);
            fs.applyBC_MBN(A, Direct.All, 3);
            fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[]{topBoundary}, null, 3);

            solver.setMatrix(A);
            solver.PCGSSOR(RHS, dx, 1.5, 1);
            double[] diff = MVO.add(dx, -1.0, x);
            if(MVO.L2Norm(diff) <= 1.0e-8 * MVO.L2Norm(x)) {
                break;
            }
            Arrays.setAll(x, i -> dx[i]);

            A.clear();
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("SolExample662.dat", x);
        A.clear();
        fs.assembleStiff(new double[]{k0}, BVPType.COMMON, A);
        Arrays.fill(RHS, 0.0);
        fs.applyBC_MBN(A, RHS, Direct.All, 500, 4);
        fs.applyBC_MBN(A, RHS, Direct.All, 300, 2);
        fs.applyBC_MBN(A, Direct.All, 3);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[]{topBoundary}, null, 3);
        solver.setMatrix(A);
        solver.PCGSSOR(RHS, x, 1.5, 1);
        mesh.toTecplot("LinearSolExample662.dat", x);
    }

    public static void NewtonIteration_example662() {
        Mesh2T3 mesh = new Mesh2T3().square2D(8, 4, T3Type.Left);
        mesh.scale(0.2, 0.1);
//        mesh.displayNodes();
//        mesh.displayElements();
//        mesh.displayBoundarys();

        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A0 = new SRMatrix(fs.GetNdof());
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        //double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 0.2;
        double k1 = 100;

        ScalarFunc initX = (xy, label, param) ->500 + xy[0] * (300 - 500) / 0.2;
        double[] x = fs.valueFromFunc(initX, null);
//        mesh.setValueOnBoundary(x, 500, 4);
//        mesh.setValueOnBoundary(x, 300, 2);
        ScalarFunc topBoundary = (xy, label, param) -> 500 * (1 - 10 * xy[0] * xy[0]);
        mesh.setValueOnBoundary(topBoundary, null, x, 3);

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] k0k1Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0k1Vec, i-> k0 * k1);
        double[][] coef = {k0Vec, k0k1Vec, zeroVec, zeroVec,
                k0Vec, k0k1Vec, zeroVec, zeroVec, x};
        while(k < maxIteration) {
            A0.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_1, A0);
            double[] R = MVO.add(RHS, -1.0, A0.mul(x));

            A.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_2, A);
            fs.applyBC_MBN(A, R, Direct.All, 0.0, 2, 3, 4);

            solver.setMatrix(A);
            solver.PCGSSSOR(R, dx, 1.5, 1);

            //5. 计算整体解
            Arrays.setAll(x, i -> x[i] + dx[i]);

            //6. 检验收敛性
            System.out.println(MVO.L2Norm(dx));
            System.out.println(MVO.L2Norm(R));
            System.out.println(MVO.L2Norm(RHS));
            if(MVO.L2Norm(dx) < 1.0e-10 * MVO.L2Norm(x) && MVO.L2Norm(R) < 1.0e-8) {
                break;
            }
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("FinalSolNewton.dat", x);

    }

    public static void NewtonIteration_example662_2DQ4() {
        Mesh2Q4 mesh = new Mesh2Q4().square2D(8, 4);
        mesh.scale(0.2, 0.1);
//        mesh.displayNodes();
//        mesh.displayElements();
//        mesh.displayBoundarys();

        FES2Q41 fs = new FES2Q41(mesh);
        SRMatrix A0 = new SRMatrix(fs.GetNdof());
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        //double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 0.2;
        double k1 = 100;

        ScalarFunc initX = (xy, label, param) ->500 + xy[0] * (300 - 500) / 0.2;
        double[] x = fs.valueFromFunc(initX, null);
//        mesh.setValueOnBoundary(x, 500, 4);
//        mesh.setValueOnBoundary(x, 300, 2);
        ScalarFunc topBoundary = (xy, label, param) -> 500 * (1 - 10 * xy[0] * xy[0]);
        mesh.setValueOnBoundary(topBoundary, null, x, 3);

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] k0k1Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0k1Vec, i-> k0 * k1);
        double[][] coef = {k0Vec, k0k1Vec, zeroVec, zeroVec,
                k0Vec, k0k1Vec, zeroVec, zeroVec, x};
        while(k < maxIteration) {
            A0.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_1, A0);
            double[] R = MVO.add(RHS, -1.0, A0.mul(x));

            A.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_2, A);
            fs.applyBC_MBN(A, R, Direct.All, 0.0, 2, 3, 4);

            solver.setMatrix(A);
            solver.PCGSSSOR(R, dx, 1.5, 1);

            //5. 计算整体解
            Arrays.setAll(x, i -> x[i] + dx[i]);

            //6. 检验收敛性
            System.out.println(MVO.L2Norm(dx));
            System.out.println(MVO.L2Norm(R));
            System.out.println(MVO.L2Norm(RHS));
            if(MVO.L2Norm(dx) < 1.0e-10 * MVO.L2Norm(x) && MVO.L2Norm(R) < 1.0e-8) {
                break;
            }
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("FinalSolNewton.dat", x);
    }

    public static void NewtonIteration_HeatTransfer_2DQ4() {
        Mesh2Q4 mesh = new Mesh2Q4().square2D(8, 4);
        mesh.scale(0.18, 0.1);
//        mesh.displayNodes();
//        mesh.displayElements();
//        mesh.displayBoundarys();

        FES2Q41 fs = new FES2Q41(mesh);
        SRMatrix A0 = new SRMatrix(fs.GetNdof());
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        //double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 0.2;
        double k1 = 2e-3;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] k0k1Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0k1Vec, i-> k0 * k1);

        ScalarFunc initX = (xy, label, param) ->500 + xy[0] * (300 - 500) / 0.18;
        double[] x = fs.valueFromFunc(initX, null);
        double[][] coef = {k0Vec, k0k1Vec, zeroVec, zeroVec,
                k0Vec, k0k1Vec, zeroVec, zeroVec, x};
//        mesh.setValueOnBoundary(x, 500, 4);
//        mesh.setValueOnBoundary(x, 300, 2);
        while(k < maxIteration) {
            A0.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_1, A0);
            double[] R = MVO.add(RHS, -1.0, A0.mul(x));

            A.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_2, A);
            fs.applyBC_MBN(A, R, Direct.All, 0.0, 2, 4);

            solver.setMatrix(A);
            solver.PCGSSSOR(R, dx, 1.5, 1);

            //5. 计算整体解
            Arrays.setAll(x, i -> x[i] + dx[i]);

            //6. 检验收敛性
//            System.out.println(MVO.L2Norm(dx));
//            System.out.println(MVO.L2Norm(R));
//            System.out.println(MVO.L2Norm(RHS));
            if(MVO.L2Norm(dx) < 1.0e-10 * MVO.L2Norm(x) && MVO.L2Norm(R) < 1.0e-10) {
                break;
            }
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("FinalSol.dat", x);

    }

    public static void NewtonIteration_HeatTransfer() {
        Mesh2T3 mesh = new Mesh2T3().square2D(8, 4, T3Type.Left);
        mesh.scale(0.18, 0.1);
//        mesh.displayNodes();
//        mesh.displayElements();
//        mesh.displayBoundarys();

        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A0 = new SRMatrix(fs.GetNdof());
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        //double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 0.2;
        double k1 = 2e-3;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] k0k1Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0k1Vec, i-> k0 * k1);

        ScalarFunc initX = (xy, label, param) ->500 + xy[0] * (300 - 500) / 0.18;
        double[] x = fs.valueFromFunc(initX, null);
        double[][] coef = {k0Vec, k0k1Vec, zeroVec, zeroVec,
                           k0Vec, k0k1Vec, zeroVec, zeroVec, x};
//        mesh.setValueOnBoundary(x, 500, 4);
//        mesh.setValueOnBoundary(x, 300, 2);
        while(k < maxIteration) {
            A0.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_1, A0);
            double[] R = MVO.add(RHS, -1.0, A0.mul(x));

            A.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_POLYNORMIAL_2, A);
            fs.applyBC_MBN(A, R, Direct.All, 0.0, 2, 4);

            solver.setMatrix(A);
            solver.PCGSSSOR(R, dx, 1.5, 1);

            //5. 计算整体解
            Arrays.setAll(x, i -> x[i] + dx[i]);

            //6. 检验收敛性
//            System.out.println(MVO.L2Norm(dx));
//            System.out.println(MVO.L2Norm(R));
//            System.out.println(MVO.L2Norm(RHS));
            if(MVO.L2Norm(dx) < 1.0e-10 * MVO.L2Norm(x) && MVO.L2Norm(R) < 1.0e-10) {
                break;
            }
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("FinalSol.dat", x);

    }

    public static void DirectIteration_membrane() {
        Mesh2T3 mesh = new Mesh2T3().square2D(8, 8, T3Type.Left);
        mesh.scale(0.5, 0.5);
//        mesh.displayNodes();
//        mesh.displayElements();
//        mesh.displayBoundarys();

        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 1.0;
        double k1 = -0.4;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] dxVec = new double[fs.GetNdof()];
        Arrays.setAll(dxVec, i-> k1);
        double[][] coef = {k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec,
                           k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec, x};
        fs.assembleSource(new double[]{1.0}, BVPType.COMMON, RHS);
        while(k < maxIteration) {
            fs.assembleStiff(coef, BVPType.NONLINEAR_DERIVATIVE_1, A);
            fs.applyBC_MBN(A, RHS, Direct.All, 0, 2, 3);

            solver.setMatrix(A);
            solver.PCGSSOR(RHS, dx, 1.5, 1);
            double[] diff = MVO.add(dx, -1.0, x);
            if(MVO.L2Norm(diff) <= 1.0e-8 * MVO.L2Norm(x)) {
                break;
            }
            Arrays.setAll(x, i -> dx[i]);

            A.clear();
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("membraneSol.dat", x);
    }

    public static void DirectIteration_membrane_2DQ4() {
        Mesh2Q4 mesh = new Mesh2Q4().square2D(32, 32);
        mesh.scale(0.5, 0.5);
//        mesh.displayNodes();
//        mesh.displayElements();
//        mesh.displayBoundarys();

        FES2Q41 fs = new FES2Q41(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 1.0;
        double k1 = 0.2;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] dxVec = new double[fs.GetNdof()];
        Arrays.setAll(dxVec, i-> k1);
        double[][] coef = {k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec,
                k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec, x};
        fs.assembleSource(new double[]{1.0}, BVPType.COMMON, RHS);
        while(k < maxIteration) {
            fs.assembleStiff(coef, BVPType.NONLINEAR_DERIVATIVE_1, A);
            fs.applyBC_MBN(A, RHS, Direct.All, 0, 2, 3);

            solver.setMatrix(A);
            solver.PCGSSOR(RHS, dx, 1.5, 1);
            double[] diff = MVO.add(dx, -1.0, x);
            if(MVO.L2Norm(diff) <= 1.0e-8 * MVO.L2Norm(x)) {
                break;
            }
            Arrays.setAll(x, i -> dx[i]);

            A.clear();
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("membraneSol.dat", x);
    }

    public static void DirectIteration_membrane_2DQ9() {
        Mesh2Q9 mesh = new Mesh2Q9().square2D(32, 32);
        mesh.scale(0.5, 0.5);
        mesh.displayNodes();
        mesh.displayElements();
        mesh.displayBoundarys();

        FES2Q91 fs = new FES2Q91(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 1.0;
        double k1 = 0.2;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] dxVec = new double[fs.GetNdof()];
        Arrays.setAll(dxVec, i-> k1);
        double[][] coef = {k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec,
                k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec, x};
        fs.assembleHeatSource(new double[]{1.0}, BVPType.COMMON, RHS);
        while(k < maxIteration) {
            fs.assembleHeatStiff(coef, BVPType.NONLINEAR_DERIVATIVE_1, A);
            fs.applyBC_MBN(A, RHS, Direct.All, 0, 2, 3);

            solver.setMatrix(A);
            solver.PCGSSOR(RHS, dx, 1.5, 1);
            double[] diff = MVO.add(dx, -1.0, x);
            if(MVO.L2Norm(diff) <= 1.0e-8 * MVO.L2Norm(x)) {
                break;
            }
            Arrays.setAll(x, i -> dx[i]);

            A.clear();
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("membraneSolQ9.dat", x);
    }

    public static void NewtonIteration_Membrane() {
        Mesh2T3 mesh = new Mesh2T3().square2D(32, 32, T3Type.Left);
        mesh.scale(0.5, 0.5);
//        mesh.displayNodes();
//        mesh.displayElements();
//        mesh.displayBoundarys();

        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A0 = new SRMatrix(fs.GetNdof());
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 1.0;
        double k1 = -0.2;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] dxVec = new double[fs.GetNdof()];
        Arrays.setAll(dxVec, i-> k1);
        double[][] coef = {k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec,
                k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec, x};
        fs.assembleSource(new double[]{1.0}, BVPType.COMMON, RHS);
        while(k < maxIteration) {
            A0.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_DERIVATIVE_1, A0);
            double[] R = MVO.add(RHS, -1.0, A0.mul(x));

            A.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_DERIVATIVE_2, A);
            fs.applyBC_MBN(A, R, Direct.All, 0.0, 2, 3);

            solver.setMatrix(A);
            solver.PCGSSSOR(R, dx, 1.5, 1);

            //5. 计算整体解
            Arrays.setAll(x, i -> x[i] + dx[i]);

            //6. 检验收敛性
//            System.out.println(MVO.L2Norm(dx));
//            System.out.println(MVO.L2Norm(R));
//            System.out.println(MVO.L2Norm(RHS));
            if(MVO.L2Norm(dx) < 1.0e-10 * MVO.L2Norm(x) && MVO.L2Norm(R) < 1.0e-10 * (1 + MVO.L2Norm(RHS))) {
                break;
            }
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("NewtonMembraneSol.dat", x);

    }

    public static void NewtonIteration_Membrane_2DQ4() {
        Mesh2Q4 mesh = new Mesh2Q4().square2D(128, 128);
        mesh.scale(0.5, 0.5);
//        mesh.displayNodes();
//        mesh.displayElements();
//        mesh.displayBoundarys();

        FES2Q41 fs = new FES2Q41(mesh);
        SRMatrix A0 = new SRMatrix(fs.GetNdof());
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 1.0;
        double k1 = 0.2;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] dxVec = new double[fs.GetNdof()];
        Arrays.setAll(dxVec, i-> k1);
        double[][] coef = {k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec,
                k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec, x};
        fs.assembleSource(new double[]{1.0}, BVPType.COMMON, RHS);
        while(k < maxIteration) {
            A0.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_DERIVATIVE_1, A0);
            double[] R = MVO.add(RHS, -1.0, A0.mul(x));

            A.clear();
            fs.assembleStiff(coef, BVPType.NONLINEAR_DERIVATIVE_2, A);
            fs.applyBC_MBN(A, R, Direct.All, 0.0, 2, 3);

            solver.setMatrix(A);
            solver.PCGSSSOR(R, dx, 1.5, 1);

            //5. 计算整体解
            Arrays.setAll(x, i -> x[i] + dx[i]);

            //6. 检验收敛性
            System.out.println(MVO.L2Norm(dx));
            System.out.println(MVO.L2Norm(R));
            System.out.println(MVO.L2Norm(RHS));
            if(MVO.L2Norm(dx) < 1.0e-8 * MVO.L2Norm(x) && MVO.L2Norm(R) < 1.0e-8) {
                break;
            }
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("NewtonMembraneSol.dat", x);

    }

    public static void NewtonIteration_Membrane_2DQ9() {
        Mesh2Q9 mesh = new Mesh2Q9().square2D(16, 16);
        mesh.scale(0.5, 0.5);
//        mesh.displayNodes();
//        mesh.displayElements();
//        mesh.displayBoundarys();

        FES2Q91 fs = new FES2Q91(mesh);
        SRMatrix A0 = new SRMatrix(fs.GetNdof());
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        double[] x = new double[fs.GetNdof()];

        double[] dx = new double[fs.GetNdof()];
        NewIterSSolver solver = new NewIterSSolver();
        int maxIteration = 100000;
        int k = 0;

        double[] zeroVec = new double[fs.GetNdof()];

        double k0 = 1.0;
        double k1 = 0.2;

        double[] k0Vec = new double[fs.GetNdof()];
        Arrays.setAll(k0Vec, i-> k0);
        double[] dxVec = new double[fs.GetNdof()];
        Arrays.setAll(dxVec, i-> k1);
        double[][] coef = {k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec,
                k0Vec, dxVec, zeroVec, zeroVec, dxVec, zeroVec, zeroVec, x};
        fs.assembleHeatSource(new double[]{1.0}, BVPType.COMMON, RHS);
        while(k < maxIteration) {
            A0.clear();
            fs.assembleHeatStiff(coef, BVPType.NONLINEAR_DERIVATIVE_1, A0);
            double[] R = MVO.add(RHS, -1.0, A0.mul(x));

            A.clear();
            fs.assembleHeatStiff(coef, BVPType.NONLINEAR_DERIVATIVE_2, A);
            fs.applyBC_MBN(A, R, Direct.All, 0.0, 2, 3);

            solver.setMatrix(A);
            solver.PCGSSSOR(R, dx, 1.5, 1);

            //5. 计算整体解
            Arrays.setAll(x, i -> x[i] + dx[i]);

            //6. 检验收敛性
            System.out.println(MVO.L2Norm(dx));
            System.out.println(MVO.L2Norm(R));
            System.out.println(MVO.L2Norm(RHS));
            if(MVO.L2Norm(dx) < 1.0e-8 * MVO.L2Norm(x) && MVO.L2Norm(R) < 1.0e-8) {
                break;
            }
            k++;
        }
        System.out.println(Arrays.toString(x));
        mesh.toTecplot("NewtonMembraneSol.dat", x);

    }

}
